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We study in this work the transport properties of an impurity immersed in a granular gas under 
stationary nonlinear Couette flow. The starting point is a kinetic model for low-density granular 
mixtures recently proposed by the authors [Vega Reyes F et al. 2007 Phys. Rev. E 75 061306]. Two 
routes have been considered. First, a hydrodynamic or normal solution is found by exploiting a 
formal mapping between the kinetic equations for the gas particles and for the impurity. We show 
that the transport properties of the impurity are characterized by the ratio between the temperatures 
of the impurity and gas particles and by flve generalized transport coefficients: three related to the 
momentum flux (a nonlinear shear viscosity and two normal stress differences) and two related to 
the heat flux (a nonlinear thermal conductivity and a cross coefficient measuring a component of 
the heat flux orthogonal to the thermal gradient). Second, by means of a Monte Carlo simulation 
method we numerically solve the kinetic equations and show that our hydrodynamic solution is valid 
in the bulk of the fluid when realistic boundary conditions are used. Furthermore, the hydrodynamic 
solution applies to arbitrarily (inside the continuum regime) large values of the shear rate, of the 
inelasticity, and of the rest of parameters of the system. Preliminary simulation results of the true 
Boltzmann description show the reliability of the nonlinear hydrodynamic solution of the kinetic 
model. This shows again the validity of a hydrodynamic description for granular flows, even under 
extreme conditions, beyond the Navier-Stokes domain. 
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I. INTRODUCTION 

The understanding of transport processes occurring in granular mixtures is still challenging. In the low- and 
moderate-density regimes the Boltzmann and Enskog equations, suitably adapted to account for inelastic collisions, 
have proven to provide an adequate framework for the study of granular flows P, 01- In particular, if the spatial 
gradients present in the system are weak, the Navier-Stokes (NS) constitutive equations for the fluxes of mass, 
momentum, and energy have been derived (with explicit expressions for the transport coefficients) for the model of 
inelastic hard spheres characterized by constant coefficients of normal restitution . Most of the early derivations 
were restricted to the quasielastic limit {aij « 1), thus assuming an expansion around Maxwellians at the same 
temperature 0, 0, H, 0, 0, Hi • However, the nonequipartition of energy becomes significant beyond the quasi-elastic 
limit, as confirmed by kinetic theory 

[13, [m, [13 , computer simulations [O, [H, [llJ3JlM^^ [S, [H , and 

real experiments (2ll. [23| . A more realistic derivation of the NS transport coefficients 2^25, 26, 27l| requires taking 
into account the nonequipartition of energy and applies for finite dissipation. The accuracy of this latter approach has 
been confirmed by computer simulations in the cases of the diffusion [2^, [2^ and shear viscosity [13, IMj coefficients. 
On the other hand, the practical applicability of the NS equations is limited to small spatial gradients, while many 
steady granular ffows do not fulfill in general this condition, due to the coupling between inelasticity and gradients 

SI'S- 

The physical situation we study in this work corresponds to a gas of inelastic hard spheres enclosed between 
two parallel walls at y = zfcL/2 moving with velocities ±J7/2 along the x-axis and kept, in general, at different 
temperatures T± [H, Tsl . [ssl . [H, [13, [IH . In the base steady state the flow velocity is along the i-axis and the 
hydrodynamic fields depend on the y variable only (planar Couette flow). This macroscopic state is characterized by 
a combined momentum and heat transport described by the pressure tensor Pij (y) and the heat flux q(y), respectively. 
A sketch of the geometry of the steady planar Couette flow for the symmetric choice = T- = is given in Fig. [T] 

Since granular matter is generally present in polydisperse form, the study of the Couette flow in the case of 
a granular mixture is an interesting problem from a fundamental and practical point of view. Needless to say, the 
general problem is quite intricate since, not only the number of parameters (masses, sizes, composition, and coefficients 
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FIG. 1; (Color online) Sketch of a planar Couette gas flow. The gas is enclosed between two infinite parallel walls located at 
y = ±L/2, moving along the a;-direction with velocities ±f//2, and kept at the temperature Tu,. 

of restitution) but also the number of transport coefficients are higher than in the monocomponcnt case. As a first 
step and to gain some insight into the general problem, in this paper we consider the tracer limit case, namely a 
binary mixture where the mole fraction of one of the components (tracer species, denoted by the label 1) is much 
smaller than that of the other component (excess species, denoted by the label 2). In this tracer limit, the state 
of the excess species is unaffected by the presence of the tracer particles and so its velocity distribution function /2 
obeys a closed Boltzmann equation in the low-density regime. In addition, the mutual collisions among the tracer 
particles can be neglected versus the tracer-gas collisions, so that the tracer velocity distribution function /i obeys a 
linear (inelastic) Boltzmann-Lorentz equation. This problem is formally equivalent to that of an impurity or intruder 
immersed in a granular gas, and this will be the terminology used in this paper. Since the impurity particle is assumed 
to be mechanically different from the gas particles, the dimensionless parameters characterizing the mixture are the 
coefficients of restitution ai2 and 022, the mass ratio mi/m2, and the size ratio a\ja2- 

Unfortunately, the complexity of the nonlinear Couette fiow makes its treatment at the level of the Boltzmann 
equation practically unattainable, even in the monocomponent case. Thus, here we will consider a model kinetic 
equation recently proposed for granular mixtures [3^. In the tracer limit, this kinetic model reduces to the same 
closed kinetic equation for the excess species as considered in Ref. [l^l plus a Boltzmann-Lorentz-like kinetic equation 
for the impurity particle. The kinetic equation for /2 admits an exact solution for the steady planar Couette flow [3^ . 
Exploiting the formal analogy between the kinetic equations for /i and /2, we find in this paper an exact solution 
for /i . This solution allows us to obtain the most relevant velocity moments of /i , which are directly related to the 
momentum and heat fluxes associated with the impurity. In particular, as expected, the impurity temperature clearly 
differs from the granular temperature of the gas particles, showing again the breakdown of the energy equipartition 
in nonequilibrium states. 

The exact solution found here qualifies as a "normal" or hydrodynamic solution since /i and /2 depend on space 
only through an explicit functional dependence on the hydrodynamic fields. This hydrodynamic description applies 
even at strong dissipation (i.e., beyond the quasi-elastic limit) and strong inhomogeneity (i.e., beyond the NS domain), 
as long as the densest regions of the system remain sufficiently dilute and the molecular chaos assumption holds. This 
provides a counter-example against the speculation that a hydrodynamic description for granular flows is limited 
to weak dissipation and/or weak inhomogeneities. In order to assess the reliability of this hydrodynamic solution, 
we have also solved the model kinetic equation by means of Monte Carlo simulations with Couette-flow boundary 
conditions. Comparison with the hydrodynamic solution shows that the latter is not a mathematical artifact but 
applies in the bulk region of the system, where boundary effects are negligible. This agreement between theory and 
simulations holds for system sizes L as small as a few mean free paths. 

In order to gain some insight into the expected hydrodynamic fields in the Couette problem, let us consider first a 
monocomponent granular gas. In this case, the exact energy and momentum balance equations yield 




(1.1) 
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0, (1.2) 



= 0, (1.3) 



where d ~ 2 and 3 for hard disks and spheres, respectively, n is the number density, and C, is the coohng rate due to 
the inelastic character of collisions. By dimensional analysis in the dilute limit, C, = vC,*{a)^ where v ex nT^^^ is an 
effective collision frequency for hard spheres. Equations (|1.1[) - (|1.3[) do not constitute a closed set of equations for the 
hydrodynamic fields n(y), T{y), and Ux{y), unless the constitutive equations expressing the fluxes as functionals of 
the hydrodynamic fields are known. For illustration, let us assume for the moment that the hydrodynamic gradients 
arc small enough as to justify the use of the NS constitutive equations. Due to the geometry of the problem, at NS 
order we have P^x = Pyy = Pzz = P pol . l4l| . from which, with (|1.3p . it immediately follows that the hydrostatic 
pressure p = nT is constant, i.e., 

p = const. (1-4) 
Moreover, the NS constitutive equations imply that Qx = <lz = and 

D dT dn 

P-y--^^^ (1-5) 

where 77 = (jj / v)!]* (a) is the shear viscosity, k = (p/miy)K*{a) is the thermal conductivity (m being the mass of a 
particle), and /x = {T'^ / mv) fjL* (a) is a transport coefficient absent in the elastic case (a = 1). The explicit form of the 
dimensionless functions C*(q^)i 'n*{'^)i a-^d yU*(a) is known [13, EH- Insertion of Eqs. (|1.4p and p.5p into Eqs. 

(HHI) and dLll) yield 

1 du^ 

a — — - — = const, (1.6) 
v ay 



1 fl d^^ 



2m \ v dy 



T = --/ = const. (1.7) 



Therefore, according to the NS description, the local shear rate du^/dy scaled with the local collision frequency 
v cx pjT^I"^ is a constant, and the temperature profile is such that {y~^dyf'T is a constant that depends on the 
reduced shear rate a and the coefficient of restitution a. The set of NS base steady states in the system have been 
analytically solved in a recent work [i^l . 

As said before, due to inelasticity these steady states do not have small spatial gradients (except for a w 1 [i^l) and 
thus a kinetic description beyond the NS domain is in general required in order to properly describe granular Couette 
flows. Specifically, this is even more necessary if 7 > (and this happens when the viscous heating dominates over 
collisional cooling [l^l), since in this case the Knudsen number is always greater than the one for 7 < [i^l. Such 
a description of the Couette flow beyond the NS domain was carried out in Rcf. [STj for a monocomponent granular 
gas with 7 > 0. Interestingly enough, this solution shares with the NS description the structure of the hydrodynamic 
profiles (|1.4p . (|1.6p . and (|1.7p . However, in the constitutive equations, the transport coefficients and the parameter 7 
are nonlinear functions of the shear rate a [SJ. At the same time, the solution is also able to capture normal stress 
differences (Pxx Pw 7^ Pzz) and the component of the heat flux along the flow direction (^j, ^ 0), which are all 



nonlinear effects l37| . All theoretical results in Ref. [37[ compare well with Monte Carlo simulations of the Boltzmann 
equation, showing the reliability of the kinetic model beyond the quasi-elastic limit. As an illustrative example of the 
necessity of a nonlinear descriptio n, we briefly an alyze the case 7 = 0, which occurs for a threshold value of a that in 
the NS description is a^^^(Q!) = ^ d(^*{a)l2r\*{a) and in the nonlinear Couette flow is ath(a) = \/c?C*(a)/2[l + C*(q^)] 
[37| . We show in Fig. [5] the disagreement between both values, which becomes very apparent for values far from the 
quasielastic limit. As shown in Ref. [s^ . the nonlinear prediction at\i{a) agrees very well with Monte Carlo simulations 
of the Boltzmann equation. 

We propose in this work a theoretical solution of the nonlinear hydrodynamic profiles for the impurity that exhibits 
absence of mutual diffusion (i.e, flow velocities are equal for impurity and excess components). Furthermore, this 
solution for the impurity also obeys equations of the form (jl.4p . (|1.6p . and (|1.7p . We will use a numerical solution of 
the kinetic equation by a Monte Carlo method in order to show that the theoretical solution we propose matches the 
hydrodynamic profiles and transport coefficients that result from the kinetic equation. Furthermore, with the use of 



4 



0.8 




a 



FIG. 2: (Color online) Plot of the threshold value of the reduced shear rate, ath(") for a three-dimensional granular gas in the 
planar Couette flow. The dashed line is the result a^^(a) = \/d^^*{7^t)jQ/rfl^ obtained from the NS equations, while the solid 
line is the prediction ath(Q) = a/ d(^* (a) /2[1 + C*('^)] from an exact solution of a kinetic model of the Boltzmann equation [37| . 
The separation between both curves is a measure of the limitations of the NS description. 

the numerical solution we show that the hypotheses, notably the absence of mutual diffusion, used in order to find 
our hydrodynamic solution are actually always true in a wide range of system parameters (including shear rate and 
inelasticity). In addition, we present preliminary Monte Carlo simulations of the Boltzmann equations which confirm 
the hydrodynamic profiles predicted by the nonlinear hydrodynamic solution of the kinetic model. 

This paper is organized as follows. The kinetic model for the mixture is described in Sec. |TT1 Then, the physical 
problem we are interested in is introduced. Section IIIII presents the exact hydrodynamic solution to the kinetic 
equations for /i and /2, with explicit expressions for the heat and momentum fluxes of both species. The simulation 
method is described in Sec. IIVI and comparisons between the theoretical predictions and the simulation results is 
carried out in Sec. [V] Finally, the results arc summarized and discussed in Sec. IVII 



II. KINETIC MODEL FOR GRANULAR MIXTURES 



Let us consider a mixture composed by smooth inelastic disks {d ^ 2) or spheres {d ~ 3) of masses m.^ and diameters 
Gi , the inelasticity of collisions between a sphere of species i and a sphere of species j being characterized by a constant 
coefficient of restitution < a^- < 1. We will focus on the dilute limit, i.e., the mean free path of the particles is 
much larger than their sizes. The relevant hydrodynamic fields are the number densities n^, the flow velocity u, and 
the temperature T. They are defined in terms of moments of the velocity distribution functions fi{r,-v;t) as 



dv/,(v), (2.1) 

pu = ^ miHiU, = ^ m, / dv v/i(v), (2.2) 
i i 

i i 

where V = v — u is the peculiar velocity, n = n.i is the total number density, p ~ pi ~ rniUi is the total 
mass density, and p is the pressure. Furthermore, the second equality of Eq. (|2.2[) and the third equality of Eq. (|2.3[) 
define the flow velocity Uj and the partial kinetic temperature T.^ for each species, respectively. In addition, in the 
dilute limit the pressure tensor and the heat flux qj associated with species i are given by 



P. = m, J dvVV/,(v), = ^/ d^V^VMv 



). (2.4) 
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In the low-density regime the distribution functions fi obey a set of coupled nonlinear Boltzmann equations |24| : 

(St+v-V)/, = ^J,,[v|/,,/,], (2.5) 



where Jy[v|/i, fj] denotes the inelastic Boltzmann operator that gives the rate of change of fi due to collisions with 
particles of species j. It is given by 



Jvi^ilhJA = J d?e(?-gi2)(a-gi2) 

X [a„.V,(r, v[,t)f,{r, V^,t) - Mr, vj, t)/j(r, V2, t)] 



(2.6) 



In Eq. (|2.6p . d is the dimensionality of the system, = (at + aj) /2, S is a unit vector along the line of centers, 
G is the Heaviside step function, and gi2 = vi — V2 is the relative velocity. The primes on the velocities denote the 
initial values {v']^,V2} that lead to {vi,V2} following a binary (restituting) collision: 



v'l = vi - (1 + a,,/) (5^ • gi2)5^ 



(2.7a) 



V2 = V2 + (l + a,/) (5^ • gi2)5^, 



(2.7b) 



where /iy = m^/ (m^ + m^), so that fiij + fiji = 1. 

However, due to the mathematical complexity of the Boltzmann equation, and in order to describe general nonequi- 
librium states, it is useful to replace the Boltzmann collision operator [v|/i, fj] with a more tractable model operator 
Kij [v|/i, fj] that reproduces the coUisional transfers of mass, momentum, and energy of the true inelastic Boltzmann 
operator, namely 



idv I V I J,,[V|/,;,/,] = I V I A;,[v|/,;,/,], 



(2.8) 



Extending the model proposed by Brey et al. [43j for the monocomponcnt case and enforcing Eq. (j2.8|l in the Gaussian 
approximation, we have recently proposed the following model kinetic equation for inelastic mixtures [s^ : 



1 + 



dtf. + V . V/, = - ^ ^ [/.(V) - /.,(V)] + • [(V - U,) /,(V)] \ , 



(2.9) 



where 



47r(''-i)/2 



dT{d/2) ^ ''^ I m. 



1/2 



is a velocity- independent effective collision frequency of a particle of species i with particles of species j , 



1 + ^+ (u, -u,) 



m,T, 2d T,, 



(1-4) 



(2.10) 



is the contribution to the cooling rate of species i due to the inelastic collisions with particles of species j, and 



(2.11) 



fij{w) = rii [ ) exp 



is a reference distribution function. In the above equations. 



mi 2 
-27^ ("-"-V 



T, 



drii 



dv {^^-u,ff,=T,--^{u,-uy 



(2.12) 



(2.13) 



(2.14) 
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(u.-u,-)^ 
2d 



Tj - T, 



Ti/nii + Tj/rrij 



(2.15) 



We now specialize to the problem analyzed in this paper, namely a binary mixture where one of the species (i = 1) 
is present in tracer concentration (ni/n2 0). In this case, Eqs. (|2.2p and (|2.3p imply that u = U2 and T = T2. In 
addition, the mixture is subjected to the steady Couette flow (see Fig. [T]), so that the spatial dependence of all the 
quantities is limited to the y variable. In the tracer limit, the state of the excess component (i = 2) is not disturbed 
by the presence of the impurity and so Eq. (|2.9p for i = 2 becomes 



dh 
' dy 



- 122) + ■ [(^ ^ "2) h] 



where, according to Eqs. p.l0p ~ (|2.15p . V2, C22, and /22 are given by 

l + a22 S^(''-i)/2 



V2 



-V22-, V22 



dV{d/2) 



n2at' 




(2.16) 



(2.17) 



ah. 



-V22 



1 — Q;22 



-1^2, 



(2.18) 



/22(v) = 712 



m2 
2^2 



d/2 



exp 



m2(v - U2)^ 

2To 



Taking moments in Eq. (|2.16[) . one gets the balance equations of momentum and energy in the steady state: 

dP2,xy dP2,yy 



dy 



dy 



= 0, 



(2.19) 



(2.20) 



dq2,y , du2.x 



dy ■ dy '"'"^ 2 
Since the impurity only collides with particles of the host gas, Eq, 

dfi 



C22'T-272. 



for i ~ I reduces to 



where 



1 + Q'12 



-1^12 



(2.21) 



(2.22) 



(2.23) 



and 1^12, C12, and /12 are defined by Eqs. (P7T(I)) - ((TT5)l with T2 = T2 = T. The kinetic equations and f^T^i must 

be supplemented by appropriate boundary conditions representing the relative motion of the plates at y = ±L/2. 

The main advantage of the tracer limit is that /2 obeys a closed (inelastic) kinetic equation (the same equation as 
the monocomponent granular gas). Once solved, the moments n2{y), U2(y), and T2{y) can be inserted into Eq. (|2.22[) 
to get a closed equation for fi. Despite the simplicity of the kinetic model with respect to the original Boltzmann 
equation, the search for an exact solution to the nonlinear Couette flow problem is a formidable task. In the case 
of a monocomponent gas, an exact hydrodynamic solution was found in Ref. (37| . Of course, this solution holds for 
the kinetic equation (|2.16p of the excess component. Based on this solution, in the next section we obtain an exact 
hydrodynamic solution for the kinetic equation (j2.22p of the impurity. 



III. HYDRODYNAMIC SOLUTION BEYOND NAVIER-STOKES ORDER 



A. Excess component 



As said before, an exact solution to (|2.16p was found in Rcf. [37j . Such a solution is characterized by the following 
hydrodynamic profiles; 



P2 = n2T2 = const. 



(3.1) 
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Id , ^ 

-U2.X = a = const, (3-2) 



V2{y) dy 



2m2 



5.2 



My ) dy 



T2 = ~l{a, a.22) = const, (3.3) 



where 7(0,022) > is a dimensionless nonlinear function of the shear rate a and the coefficient of restitution Q!22. 
This quantity (henceforth called thermal curvature coefficient) characterizes the curvature of the temperature profile 
as a consequence of both the viscous heating and the coUisional cooling. The form of the profiles (|3.ip - p.3p coincides 
with the profiles (|1.4p . (|1.6p . and p.7p predicted by the NS description, except that the thermal curvature coefficient 
7 differs from its NS value and is determined consistently, as shown below. The solution to Eqs. p.2p and p.3p is 

W2,x(s) = as, T2(s) = T2(0) + es - TO27s^, (3.4) 

where the scaled variable s is defined as 

s{y)^ \\y'v2{y'), (3.5) 
Jo 

and e is an arbitrary constant that vanishes if the two wall temperatures are equal but is nonzero otherwise (T+ 7^ T_ ) 

For convenience, we refer the velocities of the particles to the Lagrangian frame moving with velocity U2.x(s). In 
this frame, Eq. (|2.16p can be rewritten as 

(1 - + yyds - aVydv^ - ^C2 V • c?v) /2(s, V) = /22(s, V), (3.6) 

where 

>* _ C22 _ 1 - «22 „^ 

C2 - — - (3.7) 

and the derivative 9, is taken at constant V = v — U2(s). Note that the dependence of the reference distribution /22 
on both s and V is explicit. Taking this into account, the hydrodynamic solution to Eq. p.6p is [s^ 

/•oo 

/2(s,V)=== / dwe-(^-*'^=*)'"e-^('"'^^*'^«^=e°"'^«^^-/22(s,e^'^>V), (3.8) 
Jo 

where 

r(u;,C2)^^(e^'^^'"-l). (3.9) 

The action of the operators e""^^"^' and e^^^a^v-^ on an arbitrary function g(s, V) is 

e-^^«^^5(s,V)=g(s-TK„V), e'^"'^"^-^ g(s, I4) = g(s, 14 + a^^K,), (3.10) 



respectively. The solution (|3.8p clearly adopts the form of a hydrodynamic or normal solution since its spatial 
dependence only occurs through a functional dependence on the hydrodynamic fields n^is')^ U2(s), and T2(s). This 
provides a neat example of the existence of normal solutions beyond the NS domain. The solution (|3.8p depends 
parametrically on the shear rate a, the coefficient of restitution a22 and the thermal curvature coefficient 7. However, 
only the two first parameters are independent since, as indicated by the notation in Eq. ()3.3p . 7 is a nonlinear function 
of a and a2i- The parameter 7(0, q;22) is determined by imposing the consistency conditions 

y Aw{\y,v''}{h-h2)={^.^M- (3.11) 

While the first two conditions are identically satisfied regardless of the value of 7, the third condition in (|3.1ip leads 
to the following implicit equation [U 

C* 2fl^ 
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Here, we have introduced the mathematical functions 



Jo 



V^eiw, y, z)e^'('"'y'^)erfc (6'(w, y, z)) - 1 



(3.13) 



dy 

= -- / dw e~(i+")'"w™0(«;, z) I [1 + 29^{w, y, z)] e^'^^'^'^^erfc {e{w, y, z)) - 20{w, y,z)j, 

(3.14) 

where erfc(x) is the complementary error function and 

1 z 



e{w,y,z) 



2V2yl-e-5- 



(3.15) 



The representation p.l2p exists only for 7 > or, equivalently, for a > ath, where, as discussed in the Introduction, 
the threshold value ath of the shear rate corresponds to 7 = 0. In this case, Fo_m(0, C2) = Fi^,n{0, Q) " ^ i^^^ the 
appendix] and so 



*\2 



«th = ^C2*(l+C2) 



(3.16) 



In the case a — ath the viscous heating is exactly balanced by coUisional cooling. This state corresponds with the 
well-known simple shear flow [if e = in p.4p ] but also to a non-uniform steady flow (for e 7^ 0) that has been reported 
recently [i^ . 

Once the parameter 7 is obtained from Eq. (|3.12p . the velocity distribution function is completely determined from 
Eq. (|3.8p . Its relevant moments provide the momentum and heat fluxes. The nonzero elements of the pressure tensor 
are given by [13] 



P2 



1 



P2,yy _ 1 



P2 1 + C2* 

P2,zz 1 



P2 



3 + ^0,0(7, C2) + a' [^^0,2(7, C2*) + 2Fi,2(7, C2)] , 
+ ^^o,o(7,C2) + 2Fi,o(7,C2), 
+ F0A1X2), 



I + C2 



(3.17) 
(3.18) 
(3.19) 



■2,xy 
P2 



+ i^o,i(7,C2) + 2Fi,i(7,C2*) 



(3.20) 



The requirement [P2.XX + P2,yy + {d — 2)P2,z^]/p2 = d is equivalent to the consistency condition (|3.12[) . Equation 
p.20[) strongly differs from Newton's shearing law [see Eq. (jl.Sp ] since the quantity enclosed by square brackets in 
Eq. (|3.20p is a highly nonlinear function of the shear rate a through the thermal curvature coefficient 7. For instance, 
at a22 = 0-8 and a = 1 the magnitude of P2,xy is about half its Newtonian value. 

Next, we consider the heat flux components Q2.X and ^72,^- The latter can be easily determined in terms of P2^xy 
making use of the energy balance equation (|2.21[) . according to which q2_y is linear in s. Consequently, one gets 



92„ 



P2 



2m2V2l 



\P 



2.xy\ 
P2 



d,A dT2 
dy ■ 



(3.21) 



where we have taken into account that dsT2 is also linear in s [see Eq. (|3.4p ]. Equation (|3.2ip can be seen as a 
generalized Fourier's law in the sense that q2.y is proportional to the thermal gradient with an effective thermal 
conductivity that is a nonlinear function of the shear rate. The evaluation of the component q2^x is much more 
involved. Multiplying both sides of Eq. p.Sp by V^Vx and integrating over velocity, one gets [13] 



q2,x 



P2 



dT2 



^a[G(7,C2)+«'i?(7,C2)] . 
TO21/2V27 oy 



(3.22) 
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where we have called 



G(2/,z) 



dw e 



-(i+icj 



w 



-X{e{w,y,z)) + Y{e{w,y,z)) 



(3.23) 



H{y,z)= / dwe-^^+i'^^'>"'w^Y{d{w,y,z)). 



(3.24) 



Here, 



(3.25) 



Y{e) 



2(1 + e^) - y^e{3 + 20^)e'^ erfc {6) 



(3.26) 



The existence of a component of the heat flux orthogonal to the thermal gradient and parallel to the flow direction 
goes beyond Fourier's law. In fact, q2^x is at least of order a(dT2/dy) and so Eq. (|3.22p can be seen as a generalized 
Burnett effect. 



Impurity particle 



Once the hydrodynamic state of the excess component has been characterized, we next want to analyze the hydro- 
dynamic state of the impurity particle. 

First, some useful information can be extracted by taking moments in Eq. (j2.22p : 



dy 



= 0, 



(3.27) 



dy 



(3.28) 



dq 



du2 



dy dy 



^ni(ri-Ti2)-^(ui2-U2)^ 



2Ci2ni 



Ti-^(ui-U2)2 



(3.29) 



Next, wc guess (to be confirmed later) that the hydrodynamic state of the impurity is enslaved to that of the granular 
gas in the sense that 

(i) there is no mutual diffusion, i.e., VLi{y) = \i2{y), 

(ii) the mole fraction ni{y)/n2{y) is uniform, and 

(iii) the temperature ratio x = Ti{y)/T2{y) is also uniform. 

The latter parameter x must be a function of the mass and size ratios 

mi fTi 



1712 



0-2 



(3.30) 



the reduced shear rate a, and the coefficients of restitution q;22 and ai2. Of course, the temperature ratio is x = 1 
when the impurity is mechanically equivalent to the gas particles (p, = cu = 1, ayi = 0^22)- Taking into account 
assumption (i), Eqs. (I3.28|) and (|3.29|) become 



dy 



0, 



(3.31) 



aPi,xv = -7:niTi — 1 - -— + - — 



1^2 



12 



,12 



Ti 



(3.32) 
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Furthermore, assumptions (ii) and (iii) imply that the product niTi and the ratios T12/T1, 1^1/1^2, and Ci2/^^i are 
constant quantities. The three latter are given by 



_ 2^.(1 - x) 

Ti (1+m)^x' 



(3.33) 



V2 



1 + Q;i2 I I + 

1 + Q!22 



d-1 



M + X 

2fi 



(3.34) 



Ci 



C12 ^ A^ + X 
(1 + Ai)'x 



(1 - ai2). 



(3.35) 



From a formal point of view, the kinetic equation (|2.16[) becomes Eq. p.22[) by making the changes /2 fi, 
/22 ^ /12, 1^2 ^ '^i and C22 Ci2- The formal change /22 /12 implies the changes n2 — > ni, TO2 — > mi, and 
T2 Ti2. It is then convenient to introduce the auxiliary quantities 



19 1^2 

vi{y)dy vi 



(3.36) 



7 



27711 



1 d 



Tr 



Ti fj. 



7- 



(3.37) 



Equations (j3.36p and (j3.37p . along with 711T12 = const, define the profiles of the fields characterizing the distribution 
hmction /12. 

The formal mapping described above allows us to easily get the moments of /i from comparison with those of /2. 
In particular, the two first self-consistency conditions are verified, namely 



j dv{l,V}(/i-/i2)-{0,0}, 



regardless of the values of 7 and x- The third self-consistency condition reads 

^/dvT/^(A-A2)=niTi(l-^ 



(3.38) 



(3.39) 



This condition determines the temperature ratio x- To evaluate the left-hand side of Eq. (|3.39p . it is convenient to 
obtain first the nonzero elements of the partial pressure tensor Pi. Taking into account Eqs. (|3.17[) - (|3.20[) . one gets 



Pi 



1 



niT] 



12 



+ 2 ^ 

1 + Ci (1 + Ci)^ 



+ i^o,o(7, Ci) + i^o,2(7, Ci) + 2Fi,2(7, Ci) 



(3.40) 



P^ y y ^ 1^ 
niTi2 l + Q-i 



+ J^o,o(7,Ci) + 2Fi,o(7,Ci), 



(3.41) 



Pl.zz _ 1 

niTi2 l + li 



+ -F^o,o(7,Ci), 



(3.42) 



(1 + Ci)^ 



a^^o,i(7,Ci)-2aFi,i(7, Ci), 



(3.43) 



where the functions -Fo.m(j/, ^) and i^i_m(y, z) are defined by Eqs. (|3.13p and (|3.14p . respectively. Condition (|3.39|) is 
equivalent to Pi^xx + Pi,yy + {d — 2)Pi^zz — dniTi, yielding 



22^ 



Ti 



T12 1 + Ci/ (i + Ci)=^ 



2i^i,o(7, Ci) + rfi^o,o(7, Ci) + a' 2^1,2(7, Ci) + ^^0,2(7, Ci) 



(3.44) 
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For given values of the reduced shear rate a and the mechanical parameters of the system (a22, CH2, lo), Eq. 

(|3.44p . complemented with Eq. (|3.12p and the relations p.33p - (|3.37p . becomes a nonlinear closed equation for the 
temperature ratio x- which must be solved numerically. In the case of mechanically equivalent particles, Eq. (|3.44p 
yields x = 1 s-iid is equivalent to Eq. (|3.12p . Insertion of this solution into Eqs. (|3.40p - (|3.43p gives the elements of 
the pressure tensor Pi. 

We consider now the heat flux associated with the impurity. According to Eq. (j3.32p . qi^y is linear in s. Since dsT2 
is also linear in s [cf. Eq. (|3.4p ]. one can write 



91. 



\P'l,xy\ 

niTi 



d vi 

2172 



Tl2 

Ti 



m 

dy ■ 



(3.45) 



To get the x-component of the heat flux, we make use again of the formal mapping described above. Thus, from Eq. 
([X^ we obtain 



qi.x 



12 



G(7,Ci)+a'i?(7,Ci) 



dTi2 
dy 



(3.46) 



C. Generalized transport coefficients for the impurity particle 



In order to characterize the momentum and heat transport associated with the impurity particle we introduce five 
generalized transport coefficients. The shear stress Pi^xy defines a (dimensionless) nonlinear shear viscosity coefficient 
■qi as 



Pl,xy = -Vi 



niT2 du2,x 



vi dy 

The anisotropy of the normal stresses can be measured through the viscometric coefficients A^i and Mi. 



Pi, XX Pl^yy 

niTi mTi 



The heat flux defines a generalized thermal conductivity coefficient Ai and a cross coefficient 01 as 

d + 2 niT2 dT2 d + 2 niT2 8X2 

9i,a = -Ai— qi,x = (Pi^ 

2 TOiJ^i oy 2 mii'i oy 



(3.47) 



(3.48) 



(3.49) 



From Eqs. p.40p - p.43p . (|3.45p . and p.46p it is possible to identify the expressions for these five generalized transport 
coefficients. They are given by 



^12 

J-1 



(1 + Ci)^ 



i^04(7,Ci) + 2^^1,1(7, Ci) 



(3.50) 



^1 



Tr 



I (1 + Ci)= 



F^Ml, Ci) + 2Fi,2(7, Ci) - 2Fi^o(7, Ci) 



(3.51) 



A/i=-2^Fi,o(7,Ci), 



(3.52) 



d + 2 Ti 7 



(3.53) 



2 fTuV 



G(7,Ci)+a'i?(7,Ci) 



d + 2\ Ti J ^ 
Their expressions in the limit a — s- Oth are explicitly given in the appendix. 



(3.54) 
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IV. MONTE CARLO SIMULATIONS 



As said before, the exact solution to the kinetic equation (|2.22p derived in Sec. lIIII defines a normal or hydrodynamic 
solution where its spatial dependence only occurs through the hydrodynamic fields (ni, 77,2, U2, and T2) and their 
gradients. This solution is free from boundary-layer effects and formally corresponds to idealized boundary conditions 
of infinitely cold walls {T^^ 0). For more details the reader is referred to Appendix B of Ref. [13] • The important 
point is whether or not this exact solution actually describes the steady state reached by the system, in the bulk 
domain, when subject to realistic boundary conditions and for arbitrary initial conditions. To confirm this expectation, 
one needs to solve numerically the set of time-dependent kinetic equations 

^tf^ + Vy^^-u,{f,-f,2) + ^-^■iv-u^)f, ^ = l,2. (4.1) 

These equations are solved with boundary conditions at y = ±L/2 compatible with the wall values ±J7/2 and T^, and 
star ting from an arbitrary initial condition. Specifically, we have considered Maxwellian diffuse boundary conditions 
[3^ . l45j and an initial distribution of total equilibrium. The latter choice does not imply a loss of generality in the 
base steady states that are achieved in the system and only affects the transient evolution. Both species are let to 
simultaneously evolve from the initial state. It is also to be noticed that in the numerical solution of Eq. (|4.ip there 
is no a priori assumption of equal flow velocities for the two components, i.e., eventual steady-state solutions with 
ui 7^ U2 arc let to occur. However, as we will show, this actually never happens and all the steady states found are 
consistent with Ui = U2 (absence of diffusion). 

In this paper we have employed a direct simulation Monte Carlo (DSMC) method [i^, \^ to numerically solve the 
kinetic equations (|4.ip in the three-dimensional case. The DSMC method has been extensively used to solve kinetic 
equations like the Boltzmann and BGK equations and it has proven to accurately describe transport phenomena in 
elastic gases and has also successfully been extended to flows in granular gases. In the DSMC method two steps are 
taken every time interval St: the free streaming step, during which a particle with velocity v is drifted by vdt and the 
boundary conditions are applied to those particles leaving the system, and the collision step, in which ViSt collision 
pairs are randomly selected among neighbor particles, i^i being the characteristic collision frequency in the kinetic 
equation. Our method differs from the elastic case in the addition, in the free streaming step, of the drag term which 
mimics the inelasticity in the collisions. 

The distributions fi are represented by A/i particles with velocities {vfe} and positions {yk}, k = 1, . . . , A/i. The 
system is split into A4 layers / = 1, . . .Ai of width Sy = L/A4. The particles with positions belonging in layer / 
define the densities Uij, the flow velocities u^j, and the temperatures Tij of that layer. From those quantities one 
can evaluate uij and Ci2j- The free streaming and the collision steps are briefly described below. 



A. Free streaming 

In the free streaming step the positions and velocities for both components are updated with the following rules: 

yk yk+ Vk^ySt, 

Vk ^ u,j + e-i'^''^'/^Wk-uu), (4.2) 

where / is the layer the particle k belongs in. The spatial and velocity updates (|4.2p are valid as long as the particle 
does not leave the system, i.e., \yk +Vk^ySt\ < L/2. Otherwise, the particle is reentered by applying thermal boundary 
conditions. If the particle "crosses" a wall, then 

Vfc^±([//2)X + Wfc, (4.3) 

where the velocity components Wk^x, Wk,z are randomly picked from Maxwell distribution functions (at a temperature 
T^) whereas Wk^y = (upper and lower signs for top and bottom wall collision, respectively) with w > being a 
random velocity sampled from the Raylcigh probability distribution 

F,(^,) = ^e-™-^/2^". (4.4) 

The new position after wall collision is 

2/. - ±L/2 + wk.y (5t ^Id2L^\ . (4.5) 

V Vk,y j 
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B. Collision step 



For each layer / a number UijSt of particles is randomly selected among those belonging in the layer. Then the 
velocity Vj. of each one of those particles is replaced by 

Vfc^u,,/ + Vfc, (4.6) 

where is a random velocity sampled from a Maxwell probability distribution, with temperatures T12 and T2 for 
species i = 1 and i = 2, respectively. 



C. Time and length scales and simulation technical facts 



In the simulations, the quantities are reduced using £ and r as length and time units, respectively, where 



4(1 + 0122) y/2Trn2(T2 ' ' 

7T2 and vo = \J2T-ujlm2 being the average density of the gas particles and a reference thermal velocity, respectively. 

Since the aim of DSMC simulations is to solve a kinetic equation, it must be able to describe the dynamical processes 
occurring in the system at a microscopic level (46j . This means that the width layer 8y must be small compared to 
the typical microscopic length scale, determined by the mean free path ii. Similarly, the time step 8t must be small 
compared to the inverse of the collision frequency, . Also, for obtaining an ergodic simulation, the number of 
simulated particles Ni must be sufficiently large. We therefore have performed simulations, for both species, with 
Ni = 2 X 10^ particles, Sy = 2 x 10~^£, and St — 3 x 10~^t. In order to probe a nonlinear Couette flow with 7 > 
(a > flth), we have taken a wall velocity difference U ~ Wvq and a system size typically in the range L « 2-20£, which 
produces sufficiently large values of a. 

Taking into account that the relation between microscopic over hydrodynamic scales is given by the Knudsen 
number Kn, the bin Syh we pick for measurements of the hydrodynamic profiles, including transport coefficients, is 
of the order of Syh = 0.2Kn~^£i (note that, in our system, the reduced local shear rate a is the reference measure for 
the Knudsen number). This means that the measurements of the hydrodynamic properties are performed over sets 
of microscopic cells, i.e., an average over microscopic cells is taken for each set (macroscopic cell). In this way, the 
fluctuations of macroscopic magnitudes, typical in DSMC simulations, are greatly reduced and proflles are smoothed 
with no loss of resolution at a hydrodynamic scale. Prior to averaging over sets of cells, the hydrodynamic quantities 
and fluxes are obtained for each cell, by using the expressions that may be found in Ref. [1^. 

As already explained in Sees. |T] and IIIH the reduced shear rate a and the thermal curvature coefficient 7 are 
fundamental quantities in the problem. We measured these quantities by fitting the velocity and temperature profiles 
from the simulations to fourth-degree polynomials and extracting from these fits the derivatives appearing in the 
expressions (|3.2p and (|3.3p . 

An important point in DSMC simulations is the quality of the random number generator. We used for this purpose 
random generators from Intel MKL 9.1 [48l |. whose performance has been rigorously examined in technical tests. The 
DSMC code was written in C language and compiled with Intel C++ 10.0 compiler and run in 64-bit Linux machines. 



V. RESULTS 



A. Enslaving of the impurity 



The DSMC simulations described in the preceding section show that the steady state reached by the system is in 
agreement with the bulk profiles assumed in the hydrodynamic solution worked out in Sec. IIIH i.e., the pressure p2, 
the local shear rate a, and the local thermal curvature 7 are practically uniform. Moreover, the impurity properties 
are enslaved to those of the gas particles, namely the system evolves to Ui = U2, ni/n2 = const, and T1/T2 = const, 
in agreement with assumptions (i)-(iii) listed below Eq. (|3.29p . As an illustration. Fig. [3] shows simulation data of 
the pressure and velocity profiles for both the impurity and the gas particles at t « t and t w IO'^t. 
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FIG. 3: (Color online) Pressure and flow velocity profiles for the impurity (triangles) and the gas particles (circles) at f ~ r 
(open symbols) and t « lO^r (filled symbols). The system corresponds to ai2 — 022 = 0.9, m\/m2 — 2, ai/cr2 = 1, L = 3.251, 
and U — IOvq. At short times, the hydrostatic pressures pi are not constant and the flow velocities Ui are not equal, but the 
simulation quickly evolves to pi — const, p2 = const, and ui = U2, just like in the theoretical solution. We observed analogous 
evolutions in all simulations we performed, for a wide range of parameter values. 



B. Thermal curvature coefficient 



In the remainder of this section we compare the theoretical results derived in Sec. IIIII for d = 3 with the data 
obtained from our DSMC simulations. Before considering properties associated with the impurity, we first compare 
the shear rate dependence of the thermal curvature coefficient 7. Figure H] displays 7 versus for three values 
of the coefficient of restitution a22- 0^22 = 1 (elastic case), a22 = 0.9 (moderately inelastic case), and 022 = 0.8 
(quite inelastic case). It is observed that the theory compares well with the simulation results for the three values 
of a22 considered, even for strongly sheared gases. This confirms the reliability of a (non-Newtonian) hydrodynamic 
description for granular gases in the bulk domain and beyond the quasi-elastic limit, at least within the framework of 
the model kinetic equations used. It is apparent from Fig.[4]that, at a given value of the reduced shear rate a, the value 
of 7 decreases with increasing dissipation. This can be qualitatively understood by the tendency of the coUisional 
cooling to produce a concave temperature profile, while the viscous heating tends to produce a convex profile. In 
fact, both tendencies cancel each other at the threshold shear rate ath, where 7 = 6. The corresponding values for 
Q!22 ~ 0.9 and Q!22 = 0.8, are oth — 0.29 and oth ~ 0.43, respectively. As noted above, our analytical solution is not 
mathematically well defined for negative values of 7 (i.e., a < Oth, shaded region in Fig. (2). This restriction obviously 
does not apply to the simulations, which can reach states with 7 < 0. These states also include those in the absence 
of shearing (a = 0). States with a = are interesting and some cases have been studied, in the framework of the NS 
description and/or in the quasi-elastic limit [49j . 



C. Temperature ratio 



Let us study now the main properties characterizing the hydrodynamic state of the impurity. The parameter 
space of the problem is made of four (dimensionless) material quantities (the mass ratio fi = mi/m2, the size ratio 
Lu = Ox 1 02, and the coefficients of restitution ai2 and 022) plus the reduced shear rate a. For the sake of illustration, 
we will assume a common coefficient of restitution ayi = 0^22 = ol and a common size (lo = 1), so that the parameter 
space becomes three-dimensional. Furthermore, we focus on three values of /i (/i = 2, /i = 1, and ^ = 0.5) and three 
values of a (a = 1, a = 0.9, and a ~ 0.8), so that we consider nine different systems. For each one, we analyze the 




FIG. 4: (Color online) Shear rate dependence of the parameter 7 measuring the curvature of the temperature profile [see Eq. 
(|3.3p ] for 022 = 1 (solid line and circles), 022 = 0.9 (dashed line and squares), and 022 = 0.8 (dotted line and triangles). The 
symbols represent the simulation results, while the lines are the theoretical predictions given by Eq. (|3.12|l . 
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FIG. 5: (Color online) Shear rate dependence of the temperature ratio x = T1/T2 in the case of an impurity particle with 
(jj = cri/(T2 = 1, Qii = 022 = a, and = m\/m2 = 2 (dashed lines and circles) and /x = m\jmi — 1/2 (dotted lines and 
squares). The symbols represent the simulation results, while the lines are the theoretical predictions given by Eq. (|3.44|l . The 
top, middle, and bottom panels correspond to a = 1, a = 0.9, and a = 0.8, respectively. The dotted vertical lines indicate the 
location of the threshold value Oth(o). 
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FIG. 6: (Color online) Shear rate dependence of the nonlinear shear viscosity coefficient 771 [see Eq. (|3.47|l ] associated with an 
impurity particle with ui = ai/a2 = 1, an = 0^22 = ct, and /i = m\/m2 = 2 (dashed lines and circles), ^ = mi/m2 = 1 (solid 
lines and triangles), and /x = mi/m2 — 1/2 (dotted lines and squares). The symbols represent the simulation results, while the 
lines are the theoretical predictions given by Eq. p.50p . The top, middle, and bottom panels correspond to a = 1, a = 0.9, 
and a = 0.8, respectively. The dotted vertical lines indicate the location of the threshold value a^i^{a). 



dependence of the properties of the impurity on the shear rate. Note that, since u> = 1 and 012 = 022 , in the case 
/i = 1 the impurity is mechanically equivalent to the gas particles. 

First, the breakdown of energy equipartition, as measured by the temperature ratio x = T1/T2, is plotted in Fig. [5] 
A good agreement between theory and simulations is observed. The lack of energy equipartition is expected because 
of two reasons. On the one hand, the state of the system is far from equilibrium due to the shearing and thus Ti 7^ T2 
even in the elastic case (a = 1) [HO, [5l|. On the other hand, even in the homogeneous cooling state, the inelasticity 
drives the system out of equilibrium and, consequently, Ti ^ T2 [l^l- We see from Fig. [5] that the impurity has a 
higher (lower) granular temperature than the gas if it is heavier (lighter) than a gas particle. This agrees with the 
general trend observed in experiments [l^, [2^. Figure [5] also shows that, for a given value of a, the deviation of the 
temperature ratio from unity increases as the shear rate increases. Similarly, at a given value of a, the deviation x — 1 
becomes more important with increasing dissipation. 



D. Generalized transport coefficients 

Next, we explore the momentum and heat transport of the impurity, as measured by the rheological quantities 771, 
A^i, Ml, Ai, and 0i defined by Eqs. (|3.47p - (|3.49p . Figures [B]-[5] display the three transport coefficients associated with 
the pressure tensor. As in the case of x, the agreement between the theoretical predictions and the simulation results 
is very good. It is apparent that, regardless of the value of a, shear thinning effects are present, i.e., the nonlinear 
shear viscosity 771 decreases with increasing shear rate. Regarding the influence of the mass ratio, we observe that, 
for fixed values of a and a, 771 increases as the mass ratio increases. The influence of dissipation on 771 is smaller than 
that of fj,. In any case, although hardly apparent in Fig. [51 the value of 771 increases as a decreases at given fi and a. 
It is interesting to note that the ratio rji/x is practically independent of fi, although it exhibits a weak dependence 
on a. 

The viscometric coefficients A^i and Mi, which measure normal stress differences, are plotted in Figs. [7] and [51 
respectively. The shearing produces a strong anisotropy in the normal stresses: Pi^xx > n-iTi > Pi^zz > Pi,yy As 
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FIG. 7: (Color online) Shear rate dependence of the reduced normal stress difference A'^i [see Eq. (|3.48|l ] associated with an 
impurity particle with u = ai/a2 = 1, an = 022 = Q, and fi = m\/m2 = 2 (dashed lines and circles), jj, = mi/m2 — 1 (solid 
lines and triangles), and /i = mi/m2 — 1/2 (dotted lines and squares). The symbols represent the simulation results, while the 
lines are the theoretical predictions given by Eq. p.5ip . The top, middle, and bottom panels correspond to a = 1, a = 0.9, 
and Q = 0.8, respectively. The dotted vertical lines indicate the location of the threshold value 0^1^(0). 



expected, this anisotropy increases with the shear rate. While, for given a and a, the coefficient Ni increases as the 
impurity becomes heavier, the opposite happens with the coefficient Mi. With respect to the influence of a, it turns 
out that it is practically negligible in the case of A^i, while Mi decreases significantly as the system becomes more 
inelastic. 

Finally, the two transport coefficients Ai and (f>i measuring the heat flux arc plotted in Figs. [5]and[T(71 respectively. 
These coefficients are quite difficult to measure in the simulations near the threshold shear rate ath, since there the 
thermal gradient is very small. This explains the scatter of the simulation data near = a^jj. Again, theory compares 
quite well with simulations. This is rather satisfactory especially in the case of (j>i since this cross coefficient measures 
complex coupling effects between the velocity and temperature gradients, which are absent in the NS regime. Figure 
[5] shows that, analogously to what happens with rji, the generalized thermal conductivity Ai decreases with increasing 
shear rate. In contrast, the cross coefficient 0i has a non-monotonic dependence for small inelasticities. In agreement 
with the behavior found for ryi and A^i, both coefficients Ai and (pi decrease as the mass of the impurity decreases, at 
given values of a and a. As for the influence of a, the results show that Ai and (j>i increase with increasing dissipation, 
this effect being more important for a heavy impurity than for a light impurity. We have observed that the influence 
of the mass ratio on Ai and (j>i is signiflcantly inhibited when one considers the ratios Xi/x^ ^^id (pi/y^ , especially in 
the former case. A remarkable counter-intuitive feature is that the coefficient 4>i can turn out to be larger than Ai for 
sufficiently large shear rate. This effect is more notorious as the system becomes more inelastic and/or the impurity 
becomes heavier. In fact, in the cases ^ = 1 and ^ = 2 with a = 0.8, one has 4>i > Ai for any shear rate larger 
than flth. Taking into account the definitions (|3.49p . the situation (j>i > Ai implies that \qx\ > \qy\; i.e., the shearing 
induces a heat flux with a component orthogonal to the thermal gradient that is larger than the component parallel 
to the gradient. 
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FIG. 8: (Color online) Shear rate dependence of the reduced normal stress difference M\ [see Eq. (|3.48p ] associated with an 
impurity particle with lo = o"i/(72 = 1, an = 022 = a, and fj, = m\/m2 = 2 (dashed lines and circles), ^ = 7711/7712 = 1 (solid 
lines and triangles), and jj, = 7711/7712 = 1/2 (dotted lines and squares). The symbols represent the simulation results, while the 
lines are the theoretical predictions given by Eq. (|3.52p . The top, middle, and bottom panels correspond to a = 1, a = 0.9, 
and Q = 0.8, respectively. The dotted vertical lines indicate the location of the threshold value a^^^{a). 

E. Preliminary DSMC results from the true Boltzmann equation 

Thus far we have shown that the numerical solutions of the model kinetic equations (|4.1[) with realistic boundary 
conditions support the steady-state hydrodynamic solution derived in this paper for the same model beyond the 
small Knudsen number limit. However, the important question is whether or not such a generalized hydrodynamic 
description is supported by the more fundamental Boltzmann equation. Comparison between DSMC simulations of 
the Boltzmann equation and the hydrodynamic solution of the kinetic model shows that the answer is affirmative in 
the case of a monocomponent granular gas (37| . 

When an impurity particle is embedded in the host granular gas, the crucial point of the hydrodynamic solution 
worked out in section [nT] is the enslaving of the hydrodynamic fields of the impurity to those of the bath, as expressed 
by assumptions (i)-(iii) below Eq. p.29p . We have performed preliminary DSMC simulations of the Boltzmann 
equation for the host gas and the coupled Boltzmann-Lorentz equation for the impurity particle and have observed 
that the properties (i)-(iii) are indeed satisfied in the steady state and in the bulk domain. As an illustrative example, 
Fig.llllshows the pressure and velocity profiles, as obtained from DSMC simulations of the Boltzmann equation, for a 
system similar to that of Fig. [3] but with a larger separation between the plates. Again, in the steady state (and also 
practically during the transient regime) one has ui = U2. Moreover, both pi and p2 are practically constant in the 
bulk region. As in Fig. El pi/ni > P2/n2, but this effect is smaller than in Fig. [3] because now L is larger and so the 
shear rate a is smaller. Moreover, as exemplified by Figs. [3] and [Til we have observed that the boundary effects are 
more important in the case of the Boltzmann description than in that of the kinetic model. A more exhaustive study, 
including the temperature ratio and the generalized transport coefficients, is ongoing and will be published elsewhere 
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FIG. 9: (Color online) Shear rate dependence of the nonlinear thermal conductivity coefficient Ai [see Eq. (|3.49|) ] associated 
with an impurity particle with uj = axjai = 1, Qh = a22 = ct, and = mi/m2 = 2 (dashed lines and circles), /i = mi/mg = 1 
(solid lines and triangles), and ^ = m\/m2 = 1/2 (dotted lines and squares). The symbols represent the simulation results, 
while the lines are the theoretical predictions given by Eq. (|3.53p . The top, middle, and bottom panels correspond to a = 1, 
a = 0.9, and a = 0.8, respectively. The dotted vertical lines indicate the location of the threshold value atij(a). 



VI. CONCLUSIONS 



In this paper we have analyzed the transport properties of impurities immersed in a granular gas under nonlinear 
steady planar Couette flow. We have focused on situations where the shear rate is large enough as to make the 
viscous heating term prevail over the inelastic cooling term in the energy balance equation. In these conditions the 
NS description is in general inadequate, as illustrated by Fig. [21 and so a more fundamental kinetic theory is needed. 
Due to the mathematical complexity of the Boltzmann equation, here we have used a kinetic model for granular 
mixtures recently proposed by the authors [1^ . Our approach differs from a recent work [sll on a bidisperse granular 
fluid under Couette flow, where a continuum description is used. In addition, the present work extends to inelastic 
collisions a previous study [sot carried out for ordinary gaseous mixtures. 

Two different and complementary routes have been considered. First, an exact hydrodynamic (or "normal") solution 
for the steady state has been found. This solution applies for arbitrarily large shear rates a (larger than the threshold 
value Oth corresponding to the simple shear flow) and arbitrary values of the parameters of the system (coefficients 
of restitution, masses, and sizes). Progress has been made taking advantage of a formal mapping between the kinetic 
equation for the gas particles (whose exact hydrodynamic solution was found in Ref. [sj) and the kinetic equation for 
the impurity. This formal mapping is possible once it is guessed that the hydrodynamic profiles of the impurity are 
enslaved to those of the gas particles, i.e., n\jni and TxjTi are uniform and Ui = U2 (no diffusion). Second, wc have 
solved the set of two coupled kinetic equations by means of a DSMC method [i^, with realistic boundary conditions. 
The numerical solution shows, in the context of our kinetic model description, the validity of the assumptions we 
make in the calculation of the theoretical solution. Furthermore, we have not found ranges of parameter values where 
these assumptions are not accurately fulfilled in the bulk of the fluid. Thus, an important corollary of this work is 
that under Couette flow and for our kinetic model the impurity never shows steady-state diffusion with respect to the 
granular gas where it is immersed (even in a strongly sheared system). 

In order to characterize the nonequilibrium state of the impurity, we have selected a number of relevant dimcnsionless 
coefficients. The basic one is the temperature ratio x — Ti/T2, quantifying the lack of energy equipartition between 



20 



2.5 
2.0 
1.5 

0.5 



2.0 

1.5 h 
^-1.0 
0.5 

2.0 
1.5 

0.5 
0.0 



T ' 1 ' 1 ' 1 — 

a=l 



^ — \ — ^ 



H — ' — \ — ^ 



■J^ 



a=0.9 
• ~ ~ ■ 



^ — ^ 



H \ ^ \ K 



a=0.8 . 
.- ^ 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 

2 

a 



FIG. 10: (Color online) Shear rate dependence of the cross coefficient (j>i [see Eq. (|3.49|l ] associated with an impurity particle 
with uj = ai/a2 = 1, Qii = Q22 = a, and fi = m\/m2 — 2 (dashed lines and circles), ^ = m\/m2 — 1 (solid lines and 
triangles), and fi = mi/m2 = 1/2 (dotted lines and squares). The symbols represent the simulation results, while the lines are 
the theoretical predictions given by Eq. (|3.54|l . The top, middle, and bottom panels correspond to a = 1, a = 0.9, and a = 0.8, 
respectively. The dotted vertical lines indicate the location of the threshold value 0^1^(0). 
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FIG. 11: (Color online) Pressure and flow velocity profiles for the impurity (triangles) and the gas particles (circles) at i w r 
(open symbols) and t ~ lO'^r (filled symbols). The system corresponds to Q12 — 022 = 0.9, m\/m2 — 2, (Ji/(J2 = 1, L = 7.14?, 
and U — lOvo- The data have been obtained from DSMC simulations of the Boltzmann equation. 
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both species. The momentum flux defines three independent coefficients: the nonUnear shear viscosity ryi, Eq. (|3.47p . 
and the two viscometric coefficients (or normal stress differences) Ni and Mi, Eq. (|3.48p . Similarly, the heat flux 
defines the nonlinear thermal conductivity Ai and the cross coefficient (pi, Eq. (|3.49p . Notice that the coefficients 
A'^i, Ml, and 0i do not have counterparts at the NS level. In particular, the coefficient (jfi is interesting because it 
accounts for a component of the heat flux orthogonal to the thermal gradient, induced by the shearing. 

Comparison between the exact solution and the DSMC simulations shows a good agreement, thus indicating the 
existence of a hydrodynamic or normal solution, even under extreme conditions, beyond the NS regime. The results 
show that, in general, Ti is higher (lower) than T2 if mi is larger (smaller) than m2. Moreover, as expected, the 
deviation of the temperature ratio x from unity increases as the inelasticity and/or the shear rate increase. Concerning 
the generalized coefficients 771 and Ai, it is observed that both decrease as the shear rate increases, while they increase 
with increasing dissipation and mass ratio mi/m2- As expected, the anisotropy of the normal stresses increases as 
the shear rate increases. In addition, as the impurity becomes heavier, the difference between the xx and yy stresses 
increase, while the difference between the zz and yy stresses decrease. The latter effect is also present when the 
system becomes more inelastic. Finally, in general, the cross coefficient 4>i does not present a monotonic dependence 
on the shear rate. However, like in the cases of 771 and Ai, the coefficient 0i increases as the mass of the impurity 
and/or dissipation increase. Interestingly, the latter effect is so remarkable that 4>i can be larger than Ai (and hence 
\<lx\ > \<ly\) if the impurity is sufficiently massive or the system is sufficiently inelastic. 

The work carried out in this paper can be extended along several lines. On the one hand, since the states considered 
here have been restricted to conditions where 7 > (a > Oth), it would be desirable to extend the analysis to the 
complementary situations where 7 < (a < ath, shaded region in Fig. [5]). While the simulation method does not 
present any technical difficulty in the latter case, the analytical solution found in this paper involves [see, for 
instance, Eqs. (|3.12p - (|3.15p ] and so is not mathematically well defined when 7 < 0. However, we have observed 
that an analytical continuation of the solution accounts well for the simulation results for a range of negative values 
of 7 [5^. Another possible alternative to overcome this technical difficulty is to carry out a perturbation solution 
in powers of 7, exploiting the fact that I7I is a small parameter in the region a < ath, as preliminary computer 
simulation results show. A second interesting problem is the extension of the tracer limit results derived here to a 
general bidisperse mixture with arbitrary composition. The main idea would be to guess hydrodynamic profiles for 
the mixture similar to those of a monodisperse system [37| . along with a common flow velocity and uniform mole 
fractions and temperature ratios. Finally, the theoretical results predicted by the kinetic model will be confronted 
with those obtained by DSMC simulations of the true Boltzmann equation. Our preliminary results show that the 
good agreement found in the monodisperse case [s^l extends to the case of mixtures, at least at a semi-quantitative 
level. 
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APPENDIX: TRANSPORT PROPERTIES ASSOCIATED WITH THE IMPURITY AT THE 

THRESHOLD SHEAR RATE 



In this Appendix we derive the explicit expressions for the transport coefficients of the impurity along the threshold 
shear rate ath(Q!22). They are obtained by taking the limit 7 — ^ 0"*" in the corresponding expressions of Sec. IIIII A 
similar study was carried out in Ref. by applying Grad's method to the Boltzmann equation. 



First, note that when y 0^ the function 9{w, y, z) defined by E 
use of the asymptotic expansion of the complementary error function [5 



3.15P goes to infinity, so that one can make 
, i.e.. 



> 1. 



Inserting this expansion into Eq. (j3.13p and performing the integral, one obtains 

Fo,,„(2/, z) « [(1 + z)-(i+") + (1 + 2z)-(i+™) - 22+"(2 + 3z)-(i+™) 



(A.i) 



(A.2) 



Since Fi^rniy,z) = ydFa.m{y,z)/dy, it follows that Fi^m{y,z) w FQ^rniy,z) to first order in y. Furthermore, the 
functions X{9) and Y{9) defined by Eqs. (|3.25p and (|3.26p . respectively, behave as 

x(^)«l, Y{e)^^, e»i. (A.3) 
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Therefore, 



G{y,z) 



[-(l + 2z)-2 + 4(2 + 3z)-2] ^2^, 



(A.4) 
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H{y, z) « — [-(1 + 2z)-4 + 16(2 + Sz)-^] y « 1. 



(A.5) 



Since both Fo_„i(7, Ci) and i^i,m(7, Ci) go to zero when 7-^0, Eq. (|3.44p becomes 

Ti 1 



T12 



1 + Ci/ (i + Ci)3 



= 0. 



(A.6) 



This is a fourth-degree algebraic equation whose physical solution gives the temperature ratio x i^i the simple shear 
flow. Once x is known, the transport coefficients are readily obtained. The coefficients associated with the momentum 
transport are, from Eqs. (|3.50p - (|3.52p . 



Ti; 



Vi = 



X 



(1 + Ci)^ 



(A.7) 



T12 



'■th 



Ti (1 + Ci)^ 



Ml = 0. 



(A., 



The evaluation of the generalized thermal conductivity Ai at 7 = from Eq. (|3.53p is trickier than before since 
substitution of Eq. (jA.7p into (j3.53p yields an indeterminate result. This difficulty is circumvented by first eliminating 
a? between Eqs. p.44p and p.53p and replacing rji by its expression (j3.50p . The result expresses Ai in terms of the 
hmctions Fo,,„ (7,(^1) and fi,m(7, Ci)- Then, the asymptotic value (|A.2p is used and the limit 7 ^ is taken. The 
final result is 



2x^ 



2 + 7Ci + Kl 



6 12 + 42Ci+37C?. 
■+2 (2 + 7Ci+6Cf)2 



^th 



(A.9) 



The limit 7 ^ of the cross coefficient (f)i is easily obtained from Eq. (|3.54p as 



d + 2 



111 

Ti 



7Ci 



(2 + 7Ci + 6Ct 



2\2 



ath 



d + 4+ 18- 



28Ci 



(2 + 7Ci+6C?)2 



(A.IO) 



where use has been made of Eqs. (jA.4p and (jA.sp . Despite the fact that there is no heat flux in the simple shear 
flow, Eqs. (jA.9p and (jA.lOp are intrinsic transport coefficients characterizing the state of the system. Equations 
(|A.7p - (|A.10p also describe the transport properties of the Couettc flow with a temperature profile linear in s. 

Equations (jA.7p-(IA.10D . when particularized to an impurity mechanically equivalent to the particles of the host 
gas, are consistent |56| with the results reported in Appendix D of Ref. [U- 
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